#include "mex.h"
#include <math.h>

static void scalarxarray(double sca, double arr[], int N)
{
    int i;
    for(i=0;i<N;i++)
        arr[i]=sca*arr[i];

}

static void neo(double* asset,int Npath,int Dtime)
{
    int i,j;
    for (j=0;j<Dtime;j++)
            for(i=0;i<Npath;i++)
                asset[j*Npath+i]=0.0;
}

static void pricex2(double* asset, double* eta, double* xi, double *X2,int Npath, int betaK, int alphaK)
{
    int i,j,k;
    //double* asset=malloc(Npath*sizeof(double));
    
// asset=X(iter,1)*ones(Npath,1);
// for jj=1:betaK
//     asset=asset+X(iter,jj+1)*eta(:,jj);
// end
// for ii=1:alphaK
//     asset=asset+X(iter,ii*(betaK+1)+1)*xi(:,ii);
//     for jj=1:betaK
//         asset=asset+X(iter,ii*(betaK+1)+1+jj)*xi(:,ii).*eta(:,jj);
//     end
// end
    
 
        for (j=1;j<betaK+1;j++)
        {
            for(i=0;i<Npath;i++)
                asset[i]=asset[i]+X2[j]*eta[(j-1)*Npath+i];
        }
        for(i=0;i<alphaK;i++)
        {
            for(k=0;k<Npath;k++)
            {
                asset[k]=asset[k]+X2[(i+1)*(betaK+1)]*xi[i*Npath+k];
            }
            for(j=0;j<betaK;j++)
            {
                for(k=0;k<Npath;k++)
                    asset[k]=asset[k]+X2[(i+1)*(betaK+1)+j+1]*xi[i*Npath+k]*eta[j*Npath+k];
            }
        }
}


void mexFunction(int nlhs,mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
        double *asset, *eta, *xi, *X2;
        int Npath,betaK,alphaK;
        int i;
        
        Npath=mxGetM(prhs[0]);
        betaK=mxGetN(prhs[0]);
        alphaK=mxGetN(prhs[1]);
//        Dtime=mxGetM(prhs[3]);
        
        plhs[0]=mxCreateDoubleMatrix(Npath,1,mxREAL);
                
//        plhs[1]=mxCreateDoubleMatrix(Dtime,1,mxREAL);
        
        asset=mxGetPr(plhs[0]);
        
        eta=mxGetPr(prhs[0]);
        xi=mxGetPr(prhs[1]);
        X2=mxGetPr(prhs[2]);
        
        
        
        
        //neo(asset,Npath,Dtime);
        pricex2(asset,eta,xi,X2,Npath,betaK,alphaK);


}
